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ABSTRACT 
A six-level baroclinic forecast model suitable for numerical 
prognosis was devised. This model includes vertical motion due to 
large-scale orographic and eddy-stress effects. The model was 
programmed in Fortran IV language and applied to a test case based 
upon initial data fields for 1200GMT, 15 December 1967. A con- 
vergent 36-hour forecast was produced. Comparisons were made with 


the verifying analyses at 12-hourly intervals. 
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1. Introduction 

One of the most successful numerical forecasting procedures is that 
of F. G. Shuman (1967) and his co-workers at the National Meteorological 
Center. This is a very sophisticated forecast scheme involving the prim- 
itive equations of motion. From the progressive improvement of NMC fore- 
cast skill-score measures, there seems to be little question that the 
NMC approach must be regarded as the direction of ultimate progress in 
the field of objective prognosis. 

One of the possible demerits to the use of the NMC system, at least 
in its present transitional stage, is the vast amount of sophisticated 
computer-time and memory required. Secondly, it employs the vertical 


1 
QT coordinate system , where Or (troposheric Y” ) is defined as 





where 
Pr = pressure at terrain 
р = pressure at Or = 0. 043 275531 
p* = pressure at the tropopause 


(assumed to be a material surface). 


In addition to the short time-step used in each NMC prediction 
iteration, a considerable amount of time is spent in interpolation from 
standard pressure levels (e.g., Fig. 1) to the O” surfaces for the pur- 
poses of the prognosis. This is followed by a final interpolation back 
to standard pressure-surface levels. Thus, it was felt here that an at- 


tempt, using vorticity methods, to initiate a quasi-hemispheric multi-level 


this description of the Shuman (f -model is not that of his most 
recent model (1967) which includes, in addition, some stratospheric 
layers. 
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forecast system could be of real advantage to the ''small-shop' forecaster. 
In this attempt the conventional FNWC-NMC grid is used, but at least 
some of the terrain effects embodied in the Shuman scheme are retained, 
even though the vorticity equation approach is used. 
In this method one expresses the horizontal wind by means of its 


usual Helmholtz expression. 


Y = K x VV + V A (2) 


Furthermore, if neither friction nor heat sources are to be considered, 
temporarily at least, Arakawa (1962) and others have shown that the 
following pair of equations represent an energy-conserving system over 


a closed mass of air: 


v^y + J(W2)+ VX: Vf = £ 2 б 
У-РУу = У% ὦ 


The use of the linear balance equation (4) would have its usual effect 
of eliminating gravity waves, initially at least. Tf terrain-induced 
friction, in its broadest sense, is to be introduced, its effect would 
impinge upon eq. (3), effecting changes primarily upon the X field, and 
to a lesser extent upon the Y field, since the balance equation is 
considered to be in effect throughout the forecast period. The local 


time derivative form of (4), 
— 2 
VüVy TV Ó; (5) 


is to be used in deriving the omega equation in section 3. 
Because of the boundary layer model which is adopted in section 2, 


it is expected that point and/or area sources of vorticity will be 
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present within the grid; therefore, the prognostic step at any time 
requires a continuous process of lateral diffusion, within the diffusion 
limits suggested by Grimminger (1941), by means of a term of the form 
r^ 
KV J to be added to the right side of (3). Diffusion also serves 
to prevent eq. (3) from developing computational instability when 
mc | 

A | Er | Ка < 
frictional considerations are injected, provided (ах) 2. 
[Richtmyer and Morton (1967)]. Неге AT is the time increment, and 


AX is the grid-mesh size. 
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2. Inclusion of large-scale friction in the vorticity equation 
If one retains the relatively simple Arakawa vorticity equation, 
while introducing friction-induced vertical motion, the resultant 


modified equation would, in the absence of orography, assume the form: 


57%), £ ) УХ: ЧЁ = = 9S, (Ke 7x T)+ K Qu (6) 


Still considering a turbulent boundary layer over a level-earth surface, 
the second term on the right side of (6) is expressible, in terms of a 


turbulence-induced vertical motion .(l , by the identity 


PD = - 9 2,(IK* Vx U) (7) 


If the Cressman approach (1960) to the parameterization of the frictional 


stress at gridpoints is employed, one may write 
LE PER Ver Vr (8) 


where Ca is actually taken from seasonal values at gridpoints after 


Kung (1963). At the top of the frictional boundary layer (indexed as 
F=1 or 2 in Fig. 1) integration of COL (р) gives the eddy-stress vertical 


motion component lg as follows: 
Lp == КУХ (Са ИЕ е) (9) 

which, after expansion, leads to 
Ne PLOEGER TE V(CaVge)* Vee] (10) 


In eqs. (8) and (9) the Kung drag coefficients are applied in connection 
with the geostrophic wind Va; at the top of the planetary boundary 
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If we revert now to a terrain-variable earth, another obvious term 


arises when orography is introduced, namely 
ἘΠῚ 97” : 
— w = — + ° — ч 

ен - 2807: νε ΨΥ Ver * Ver (11) 
This version of OL + is then combined with [fL £ of (9) so that the 
resultant maximum terrain-induced vertical velocity to be applied at 


the boundary-layer top is 

Lio + r+ (12) 
where lug is to be defined as: 
o [vr VE, - р LC s +4 vaar) VEN) (13) 


The term Var Ë V P+ qualitatively accounts for upward 
motion on the windward side of a mountain range and the reverse effect 
on the leeward side. In the computation of this term, actual values 
of pressure at terrain have been computed from consideration of the 
Berkofsky-Bertoni station elevations (1955) employing a method by means 
of which the terrain pressure is hydrostatically computed from the 
thermal character of the isobaric layer in which the terrain is embedded 
(see Appendix A). 

For the most part, the Berkofsky-Bertoni smoothed station eleva- 
tions over North America do not exceed 9000 feet. On the Asiatic con- 
tinent, however, derived terrain pressures as low as 500 mb are en- 
countered. In this regard there is a dilemma: At what level should 

War be considered representative when terrain reaches such high 


levels? ТЕ lies between 850 and 700 mb a choice of planetary 


ЕН 
boundary layer at 700 mb would appear perfectly reasonable. If Pr 
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exceeds 850 mb by 100 mb or so, the planetary boundary layer could 
reasonably be taken at 850 mb. 

Some upper limit to the magnitude of the Var vector used in Nr 
was considered appropriate in view of the calibration of Kung's sea- 
level geostrophic drag coefficients. For example, even with Pr 500 
mb, it is reasonable to consider that the effective updraft component 
need not exceed that of the 700-mb geostrophic wind over the generally 
rough terrain involved. Also, for simultaneous use of both (Lr and 

flc in (12), it is desirable to arrive objectively at a mutually 
compatible gradient level. In view of these considerations, the follow- 
ing scheme was devised to delineate the boundary layer: if Pr is greater 
than 875 mb, Рр (the boundary layer pressure) is chosen as 850 mb; if 
Pr is equal to or less than 875 mb, Pp is chosen as 700 mb. The choice 
of the planetary boundary layer which employs gridpoint wind data at 
either 850 mb or 700 mb is one of the simplifying features of the general 
frictional model introduced here. Because of the preservation of the 
orographic term, LET. ° V, with War egual either to War 
(850) or to Va F (700), it is felt that some of the main physical 
aspects of the variable atmospheric boundary layer are being preserved 
by the model, while retaining an element of convenience in a 6-level 
forecast model. Furthermore a recent study by Graham (1968) indicates 
that Kung's drag coefficient gives consistently smaller values than those 
proposed by Cressman; so their use in (8) and (9) appears reasonable with 
these geostrophic winds at the selected tops of the atmospheric boundary 
layer. 

The significance of the specification of the top of the atmospheric 


boundary layer (A.B.L.) is that this level is conventionally associated 
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with the maximum frictional vertical motion [c.f. Haltiner et al. (1963), 
Krishnamurti (1968)] which then decays in some generally exponential 
manner on either side of this level. Thus, for p # Pp» the existence of 
a frictionally-induced vertical velocity {l (p) may be postulated, in 


accordance with (13), as 
(е) = Е (2).0. T (14) 


Above the A.B.L. an analytic function Е (р) for (À (p) is chosen 


so as to imply that 


ч (100) = Q 


Fa (Pe) =/ M 


and that F (р) is a decay factor which decreases realistically with 
decreasing pressure, The choice of F (p) was made with a careful at- 
tempt at fitting the mean vertical profiles of frictionally induced 
_С) (p) depicted by Haltiner et al. (1963) in the case Ρε equal to 
850 mb. The analytic function selected for Е (р) was 
Fa (19) * |I - d ү (16) 
Фи ы 
А 
which decreases in a closely similar manner to that shown in synoptic 
evidence by Haltiner et al. (1963), and with that alluded to by Estoque 
(1957). Values of F (p) are shown in Tables la and lb. 

For elevations below the A.B.L., there is no published data to draw 
upon. Here it was desired to apply once more the limiting boundary 
conditions similar to (15) for p = 1000 mb and p = РЬ. It is also de- 
sirable to consider a function F, (P) for multiplicative association 
with Ө with a smaller rate of variation with р/р. than is the 


case with Е (р). Such a function, somewhat analogous to (16) and subject 
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to the above constraints as (16), was chosen by limiting the bracketed 


expression in (16) to the exponent 2: 


Po (1) =|1- y | 


N / ° °° ى‎ (17) 


The values of F, (P) are also in Tables la and lb. It is evident that 

| б, КОЕ | 
the formulation of $ x in (7), which will appear on the right 
side of (6), depends only upon SE. The finite-difference derivative 
with respect to pressure will in turn depend upon the scheme adopted for 
the finite-differencing of &) in the vertical, a topic to be discussed 
in connection with the solution of the omega equation. The distinction 
between bJ and Sl is that the former represents the result of partition- 
ing the total vertical motion into a component free of frictional ef- 
fects, whereas SL is the composite vertical motion resulting from 
large-scale friction. Note that the net frictional divergence in any 
column (which may include terrain) extending from the 1000-mb level to 
the 100-mb level is designed to be zero. Also since the lateral dif- 

2, 
fusion KV E is defined to be zero on the octagonal grid boundary, 
eqs. (4) and (6) still represent an energy-preserving system of equations 


in an adiabatic atmosphere. 
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з, Derivation and solution of the diagnostic omega equation 
In order to deduce the numerical solution of (6) for Ше 
| ὃ EO) 
it is clear that one reduires values for Bw at all prognostic 
levels. To obtain these vertical derivatives the diagnostic omega 
equation must first be solved at levels k = 2 to 6. 
The omega equation is derived from a combination of the vorticity 


equation and the thermodynamic equation. 


The vorticity equation, after Arakawa (1962), is now presented as 
| OW .. γὸ-ς 
yet WN VX VO no Se = t G r K VJ (18) 


Equation (18) is differentiated with respect to pressure and multiplied 
by £ to give the form 
[А 2. 
2 à | ) . OM ‘(2 F 
2 
The lateral diffusion term KV / has been omitted, since its use 
will be reserved for (18). Use of the time-differentiated linear 
balance equation in (19) leads to 
2 > Хы 
7° pg У-У Ире т рр) ТХЛ 0) 
0^ 21. 


—M ae 
The adiabatic thermodynamic equation is given as OE 


Er + Vy ‘VQ + УХ. Уфы + О(ш+П)= О αὐ 





where Ral de 

"o Pa IP 
and g is also expressible in terms of geopotential meters οσα 
(see Appendix D). Application of the Laplacian operator, Up. to 


each term of eq. (21) gives 
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70 pr +0 50) 4] V^ (vx- VB e)» V (co«o0): 0 (22) 


Then 750 сап be eliminated from eqs. (20) and (22) by subtrac- 
tion. Hence T st is та 
ζω)” OO V LJ ( 
E = YN М5») 
2 (23) 
t. - y? PL — од. 2 
Now, use of the approximations, 


2. 9 
оту, = Ды , VES VXp« £2 JO), 
together with del operator applied to (21) leads to eq. (25): 


Рур, m - τ . [V (0 δρ) vx ὧν ό(ωνΩ))] (25) 


Hence, the final form of the diagnostic omega equation is 


vito) «n 2 IO) LI SEE TEE) 
(26) 


A жа sr τὸν rol) |. von 


Equation (26) now has terms involving partial derivatives in both 
the vertical and horizontal, requiring finite-difference schemes for 
these derivatives. Appendix B presents the scheme used for πατε ο 
ferencing in the vertical, and Appendix C presents the scheme used for 
finite-differencing in the horizontal. Since GJ appears in a vertical 
derivative in eq. (3), a solution of i) at every grid point for levels 

k equal to 2 and 6 requires placing boundary conditions on omega at 
the topmost level ( k - 7) and the bottommost level (k = 1) as well as 
at the lateral boundary of the octagonal grid (Fig. 2). These conditions 


22 


were that omega would be equal to zero at all of these boundaries. 
Similarly, the Jacobian and Laplacian terms of all other parameters in 
(26) are set to zero at the lateral boundary. Once the terms of (26) 


are finite-differenced, we obtain a Helmholtz type of equation where 
a E 


Here C, is the finite-differenced result of the right side of (26), 
and р. is presumed known. Equation (26) was solved by a 3-dimensional 


Liebmann relaxation scheme where the ( n + 1 )th iterate for any point 


is given by 


T = m) 5 τ AE 
^. ia Z+ B, 


Here A is the over-relaxation coefficient and the residual g^ at 


(28) 





the nth iteration is 
^ mt^ s n E n . nel 


In (28), (29), the superscript `n indicates the nth iterate result- 
ing from sequential Liebmann solution procedure. The convergence 
criterion imposed is that when g^ becomes less than?10? for each grid- 
point, the iteration procedure ceases. For the first iteration ( t - O) 
all omegas of the right side of (29) were set to zero, and only the 

first 2 terms on the right side of eq. (26) were used to determine од, : 
When these omegas were computed for levels 2 through 6, then X. was 
determined by use of the equation of continuity in the form 


2 X cm ο πο. (30) 


ЕП units of mb 482 
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With a first determination of Т at the 5 central levels ( k = 2, ..., 
6), the diagnostic equation could now be solved for a better Ü) and 
then a better X . This WX iterative procedure is repeated until, 
for all gridpoints, at the 5 central levels the omegas reach the con- 
vergent value. 

For times other than ё = QO, the full diagnostic equation is used 
even in the first iteration. The first guess for X is that derived 
from the previous time step. 

The above procedures apply identically for levels k = 2,...,6. 

For the 1000-mb level a somewhat different procedure is necessary to 
obtain 2 for use in connection with the continuity equation (30). 


This procedure is outlined in Appendix B. This same procedure also 


applies when obtaining 20) for the prognostic step at 1000 mb. 


+ 
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4. Solution of the prognostic equation 

With the vertical derivatives of omega and the divergent part of 
the wind now available for levels k = 1 through 6, the vorticity eq. (18) 
is now ready for solution for these same levels. With all the advective 
and divergent terms of this шаша a known at time t, and after finite- 


differencing and transposition, eq. (18) leads to the Poisson-type 


differential equation 
2 E › куйры ` ' 

V (xL GA) (31) 
where G (i, j, k) represents Һе three-dimensional forcing function. 
This equation is solvable one level at a time by assuming both that Ye 
at the octagonal boundary is zero and that all terms comprising G (i, 
j, k) are also individually zero on the boundary. For solution of (31) 
one proceeds to a Liebmann two-dimensional iterative relaxation pro- 
cedure, using successive iterates in the form 

rel 1 үл 
We = We + AR (32) 
for any one level. Here, as in the relaxation process for QJ  , the 
superscript n +1 refers to the (n + 1)th iterative value, and the n- 
superscript refers to the nth iterative value. Again A refers to the 


over-relaxation coefficient and the residual R is expressed as: 
N А 4 
R ы V Wr 6; (33) 


This iterative procedure contínues until the maximum residual at any 


iteration is less than the prescribed epsilon value, е 2575 To. 


-2 
sec `, 


Because of the simplified method of determining the gradient level, 


and to take into account discontinuities in the A.B.L. occurring between 
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adjacent gridpoints, the term Kv*7 is included in (18) as a 
smoothing, or diffusion term. Several values of K were used with vary- 
ing results. The initial value was 4.0 x 107 m С as suggested by 
Danard (1966). This resulted in extreme deepening of lows and building 
of highs. The value 4.0 x 10° i een was then attempted which gave 
the opposite effect, namely to flatten out systems excessively. Also 
lateral boundary difficulties became apparent with the larger smoothing 
coefficient. The next and final value chosen was 1.2 x {οξέος which, 
while perhaps not the optimum value, gave better verification at t = 12 
hours. In accord with convergence requirements for a parabolic-type 
differential equation, the smoothing was applied to the y field of 
the time t-1 rather than time t. 

With We ( t = 0) computed initially, v for a 30 minute time- 


increment was determined by a forward time step, represented symbolical- 


ly by 
V = ү Ж А (At) (34) 


For the next time step, also for At = 30 minutes, a centered or 


leap-frog type time step was utilized with m = 1 in (35): 
mel mel yn 
Y у + Ye (2 at) (35) 


Form “> 1 all time steps were of the leap-frog type with 2 A t= 


К 


2 hours. 

New values for the stream function were thus obtained at all levels 
up to 200 mb. These values were converted to heights by the solution of 
the balance equation. In order to obtain vertical derivatives of geo- 


potential at 200 mb, and to obtain a more accurate stability parameter 
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σ (200), it was deemed necessary to obtain a good estimate of the 
100-mb height field (see di in eg. D-6). This was accomplished by use 
of regression equations developed by Lea (1961), for details see Ap- 
pendix E, 

With forecast values of Z available for levels k = 1 to 7, further 
time steps could now be computed by the procedure described in the 


diagnostic section and this section. 


27 


λα Forecast verification results 

As a preface to the discussion of the forecast verification results, 
the labeling of the maps (figs. 5-30) will be explained. A11 1000-mb 
contours are labeled in whole meters and contours are drawn at every 
60 meters. The 500- and 300-mb contours are labeled with the thousands 
digit dropped. (Thus a height of 5580 meters is labeled 580). The con- 
tours at these levels are drawn at every 120 meters. The isolines of 
vertical velocity are drawn with an interval of 2 x w millibars per 
second with the exception of the 200-mb level where contours are drawn 
at every 0.5 x 10% interval. Labels are millibars per second x ja] 

Although forecast maps were produced for / levels, only the 1000-, 
500-, and 300-mb levels are presented in order to maintain the discus- 
sion at a reasonable level. These forecasts are for 12-hour intervals 
out to 36 hours. 

The test case chosen, 1200Z 15 Dec. 1967, was one which was recent 
and thus readily available and one which shows extensive deepening and 
movement of pressure systems at most levels. 

The initial 1000-mb analysis shows an inactive 60-meter low center 
over Southern and Baja California which in 36 hours deepened to a minus 
120-meter contour and moved to the Colorado area. The 36-hour forecast 
did not verify the deepening; however, a trough-type intrusion into this 
area has occurred in the same general area. This low system aloft had 
a near vertical closed circulation throughout the atmosphere, and had a 
movement to the northeast of 600 nautical miles, with slight deepening 
at the 300-mb level. The upper level forecast moved this system north- 


ward 150 nautical miles and filled it considerably. 
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In the North Atlantic a low system of minus 120-meters, with 
troughing to the WNW at 1000 mb at initíal time, split into two minus 
120-meter systems in 36 hours: one positioned near the original low 
and the second positioned just west of Newfoundland. The 36-hour fore- 
cast for this level shows a minus 240-meter contour just to the north- 
west of Newfoundland. The initial system had a WNW tilt with height 
and a closed 8640-meter contour centered at 52N-65W at the 300-mb level, 
which in 36 hours deepened to a 8520-meter contour at 47N-60W. For the 
forecast system at 300-mb the model verified on position but filled the 
low center to a 8/00-meter contour. At the 500-mb level the prognosis 
verified with respect to position and height. 

Initially, east of Hokkaido, Japan, there was a 1000-mb low center 
with a minus 60-meter contour. This low deepened rapidly for the first 
12 hours and then the deepening subsided. In 36 hours it was positioned 
east of the Kamchatka Peninsula with a closed minus 180-meter contour. 
The surface forecast position verified exactly, but the forecast deepened 
the system too much. The initial 300-mb map indicated that a short wave 
trough was associated with this surface system; in 36 hours it had moved 
NE an amount consistent with that of the surface system while leaving a 
stationary long wave secondary trough west of Kamchatka. Both of these 
300-mb trough features verified quite well. 

At the initial time, a surface low center (minus 120-meter con- 
tour), was located off the central coast of Norway. It traveled ESE 
to Western Russia, maintaining its initial contour height. The model 
captured this movement almost exactly but again indicated too much low- 
level deepening. At the upper levels a trough associated with this 


system moved ESE to a position over the Baltic Sea in 36 hours. This 


29 


trough movement was captured well by the model at both the 500- and 300- 
mb levels and the height of the 500-mb level also verified. 

Centered at 67N-85E at the 1000-mb level at the initial time was a 
minus 180-meter contour. This low filled to a minus 60-meter contour 
in 36 hours and was positioned at /4N-110E in 36 hours. The model missed 
this system badly by deepening it throughout the 36—hour period to a 
minus 240 meter contour positioned at 72N-90E, about 250 miles to the 
west of the actual position. This system at 500 mb was only forecast 
to deepen slightly. This slight deepening and position was verified by 
the 00Z 17 December analysis. 

Aíter 36 hours, with the exception of the low system over the 
central United States, there was generally good agreement in the 1000-mb 
low systems' positions; however, the model tested produced excessive 
deepening in some of the surface low systems by an extent of the order 
of 15 mb. At the upper levels one may easily identify comparable systems 
but again the position of systems verified better than the intensities; 
and the 500-mb systems verified better than those of the 300-mb level on 
both intensity and position. Table 2 presents the pillow and root mean 
square error for all 6 levels at 12-hourly intervals based on all in- 
ternal gridpoints. 

Figures 26 to 30 show the vertical-motion or (Y charts produced 
at t = 12 hours, namely at 00Z, 16 December 1967 for the levels 850, 
700, 500, 300, and 200 mb. It is noteworthy that maximum values seem 
to be closely aligned in the vertical with extreme maximum values oc- 


curring generally at 500 mb. However, WJ was in general closer in 


300 
magnitude to Doo than was ορ” The extreme magnitudes computed were 
near 6.0 x 1077 mb вес ^ at 500 mb, which are close to values reported 


in a case study by Krishnamurti (1968). 
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6, | Remarks and conclusions 
Due to lack of sufficient core space it was not possible to main- 

tain all necessary fields for one time step in core. The choice was made 
to maintain only the information for a triad of levels at one time. This 
necessitated a large amount of reading from and writing on direct-access 
storage devices as is evidenced by the FORTRAN program. This increased 
the running time of the program considerably. In order to limit the 
continuous running time of the program, the 36-hour forecast was pro- 
duced in 3 runs of 12 hours each. The current values of height, 0) ,X , 

WY, as well as the Y at the preceding time step were stored on disk 
and these values were read to begin the first time step of the next run. 
(For the layout of direct-access devices see figures 4a, 4b, and 4c.). 
Running time was 110 minutes for the first 12 hours and 90 minutes for 
each of the remaining 2 runs. This time could be reduced to a more 
operationally acceptable figure by changing from Fortran to machine 
language programming, conversion to one-dimensional arrays, and, in 
general, increased programming efficiency. 


As originally planned, the 100-mb level was to have been completely 


inactive; that is, the equation 
2 
үу“ у, + Ур УИ иле VF =f Se» (L2 fA) (36) 


as applied at 100 mb is subject to the following boundary conditions: 


Y= 0 
Vy Vr =O 
x ¿O 
i+ O=0 
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In order to deduce improved values of vertical derivatives of the 
geopotential at 200 mb thus improving the static stability parameter, 
the heights at the 100-mb level were altered at each time step using the 
Navy regression procedure outlined in Appendix E. This procedure was in- 
tended to provide a reasonable estimate of the mean temperature in the 
200 to 100-mb layer for use in the thermodynamic equation. Due to an 
inadvertent programming error, however, the Q 100 values were converted 
to aU 100 field which varied in time. This lapse caused the altera- 
tion of the static boundary conditions above and, through the ω equa- 
tion, it affected the values of fL I, Y) in eq. (26), especial- 
ly at the upper levels. The actual desired result of this term at 200 


mb in centered finite-difference form is 


f LJ Оу, Y) зоо o] (37) 


200 тб 


The effects on this important term of the J equation were opera- 
tive from t = O but were not discernible until diverging values of LJ 
were detected at t = 42 hours, It is planned to rectify this problem and 
report the results at a future date. 

As is depicted in Table 2, the model appears to produce too much 
deepening at low levels and an overall increase in heights above 850 mb. 
The reason for this phenomenon is as yet undetermined. The fact that 
systems do indeed maintain their identity and that flow patterns, at least 
in a qualitative sense, show much similarity to the verifying analyses 
indicates that the model has some merit and further work on its refine- 


ment is justified. 
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TABLE 1 
F(p) VALUES FOR USE IN EQUATION 14 


(a) Gradient Level Assumed as 850 mb 


Pressure level F(p) 
1000 mb 0.00000 
850 mb 1.00000 
700 mb ООВ 
500 mb 0.42534 
300 mb 0.13528 
200 mb 0.03398 
100 mb 0.00000 


(b) Gradient Level Assumed as 700 mb 


Pressure level F (p) 
1000 mb 0.00000 
850 mb 0.20763 
700 mb 1.00000 
500 mb 0.56579 
300 mb 0217995 
200 mb 0.04520 
100 mb 0.00000 
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Level Diagram 


100 mb k = 7 
200 mb k = 6 
300 mb k = 5 
700 mb k = 3 


850 
875 


1000 


mb 
mb 


mb 


Gradient Level = 700 mb 





Figure 1 


Schematic representation of the pressure surfaces 
used in the diagnostic procedures. The index 

k = 1,...,7 is identified with the constant pressure 
values shown in mb. The index F = 1, 2 denotes the 
top of the atmospheric boundary level or gradient 
level, of 850 mb, or 700 mb, the selection of which 
depends upon the station pressure Pa (isj o); 
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FLOW DIAGRAM 
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(dy 22,6) 


tel 
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Figure 3a. 
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FLOW DIAGRAM (continued) 
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Figure 3b 
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TAPE 11 RECORD LAYOUT 


Record Contents 
1. 1000-mb initial Z values 
2% 850-mb initial Z values 
3. 700-mb initial Z values 
4, 500-mb initial Z values 
5. 300-mb initial Z values 
6. 200-mb initial Z values 
1% 100-mb initial Z values 
8. Z "values 
9. Z, values 

Figure 4 
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Record 


11-15. 
16-20. 
21-25. 
26-30. 


31-35. 


DISK 12 AND 13 RECORD LAYOUT 


Contents 


1000-mb Z intermediate forecast 


1000-mb Y 


1000-mb Y 
1000-mb (2 
1000-mb X 


previous time intermediate forecast 


current 


time 


intermediate 


intermediate 


contents & sequence same 


contents 


contents 


contents 


contents 


contents 


& 


& 


& 


& 


sequence 


sequence 


sequence 


sequence 


sequence 


Figure 4b 


same 


same 


same 


same 


same 
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intermediate forecast 


forecast 


forecast 


as 


as 


as 


as 


as 


as 


Records 


Records 


Records 


Records 


Records 


Records 


1-5 


1-5 


133 


Ds 


for 850 


for 700 


for 500 


for 300 


for 200 


for 100 


mb 


mb 


mb 


mb 


mb 


mb 


TAPE 14 RECORD LAYOUT 


Record Contents 

> 1000-mb Z final forecast” 

ae 1000-mb Y previous time final forecast 

з; 1000 -mb Y current time final forecast 

4. 1000-mb (Y final forecast 

5. 1000 -mb X final forecast 
6-10. Contents € sequence same as Records 1-5 for 850 mb 
11-15. Contents & seduence same as Records 1-5 for 700 mb 
16-20. Contents € sequence same as Records 1-5 for 500 mb 
21-25; Contents € sequence same as Records 1-5 for 300 mb 
26-30. Contents & sequence same as Records 1-5 for 200 mb 
31-35. Contents € sequence same as Records 1-5 for 100 mb 
36-70. Contents & sequence same as Records 1-35 for 2nd 


forecast interval 


Figure 4c 


2 

"Note: Although these records are designated as final forecast, 
they are used as input for the next computer interval run if the total 
forecast is partitioned into forecast interval segments. 
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APPENDIX A 
Computation of Terrain Pressure 


In exhibiting the method of determination of terrain pressure Pr 
the example used will be the case of terrain lower than the 850-mb 
standard height. Here 2. is the terrain height, 21 is the 1000-mb 


height and Z, is the 850-mb height, We start with the hypsometric 


2 


relationship: 
2.-2,- F4 d 
T ! 9 dis Er (A-1) 
7 
In (A-1), T is the mean temperature in the layer. Equation A-1 may also 


be expressed as 
AE > OZ Li+ 5,1 (А-2) 


where Sy is the mean specific virtual temperature anomaly as dis- 
cussed in Haltiner and Martin (1957) and subscript ST refers to stand- 
ard atmosphere value. After replacing temperature, 
WSs = + od 
Ка 2229927 — 


the following approximation for 1 + Sv is obtained in finite-differ- 


ence form: 


у = È- $, = 
(dag de 


Since À Z is known at any time, As ES 2) may be deter- 


(A-3) 


mined by the substitution of (A-3) into (A-2); hence, a value of (Zr)sr 
may be obtained. This value is compared with standard atmosphere heights 
from U, S. Standard Atmosphere (1962). When the height values are matched, 


the corrected terrain pressure is determined. 
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For cases in which the terrain extends above the 850-mb standard 
height, the same procedure is followed and sub | refers to 850 mb and 


sub 2 refers to 700 mb, and similarly for higher terrain heights. 
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APPENDIX B 
Treatment of Vertical Derivatives 


For fitting ο (vertical motion) profiles to three successive 
pressure levels (lower level subscript 1, middle level subscript 2, 
upper level subscript 3) a Lagrangian quadratic interpolation formula 


is used on omega. 


wp) = u, REED. + WEEE), EEE 1) 
(e "Pa (AR) (4 (4-3) Га NB?) 


After differentiating the above and setting p = Р. the following equa- 


tion is obtained: 
L = L), me Lucent), ως %-% (в-2) 
ФА, RER) ERNER) ORNAR) 


Similarly, differentiating (B-1) twice with respect to pressure and set- 


ting p = Р, the following equation is obtained: 


2o 2] LL LLL (472. y 23 AE S 
29 Ур, β (FR) (4 0-8) (R- P)(R- οἱ (6-%)(6-%) 


() , the frictional and terrain induced vertical motion, is treated 





the same as ÛJ in relation to the format of vertical derivatives. 

For fitting $ (geopotential) and Y (streamfunction) to three 
successive pressure levels the same procedure, as above is followed with 
the exception that ф апа Y are considered to have a quadratic 
vertical variation on a logarithmic pressure scale leading to the follow- 


ing results for d centered at p = Po: 


το 


_ + (P-D)O- P (p-P)(P-Ps) _ 4 (P-P)(P-P2) 
м") MEE” EDE №) 5: {ΡΟΝ 


where P = Ln(p) and 





оф SEE С вЫ), 3:4, (7) B-5) 
Bac EG ET 
For ° at the 1000-mb level, omega is expressed as a cubic equation 


in terms of omega at the 1000-mb (which is equal to zero), 850-mb, 700- 
mb, and 500-mb levels. After differentiating with respect to pressure 


and setting p = P1000 the following equation is obtained: 


& 
20) = +A - Po (fir fs )+ Po Ps oce + 
A (Р-Р) (6-5) (т-ру w 
= P+ +P P L) 700 + (8-6) 
(R R) (R-P) Py) 
р, P+ Pe +P, r WY evo 


A- 2 vr. P) (P3- Pa) 


Setting in appropriate values of Pg: Pi Р,, апа ZU the above becomes, 


OW — 0190476 Wt 0083333 We (B-7) 


ЫН ` — 0012857 223 


The coefficients derived above are constant with respect to time. 
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The values of the coefficients appearing in PC) and š $ 
ye 
given by eqs. (B-4), (B-5) have been stored and employed in subroutine 


for generating 


ф (Ph) k 


эр. 
δ») φῇ, 


In the uppermost level, i.e. for both $ > and Dp 6 it was found 


το... 7 


κ 


2% 9499999 99 99 6 


to be useful to make use of the Navy Weather Research Facility extra- 


polation procedure to obtain a hydrostatically-consistent best estimate 


for б 100: 
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APPENDIX C 
Treatment of Horizontal Derivatives 


The horizontal derivatives this paper is concerned with are the 
Laplacian operator, the Jacobian operator and the del-dot-del opera- 
tion. The finite difference forms of these are exhibited below; the 


superscripts refer to the designated points exhibited in Fig. 31. 


2 m^ А. + А 
Laplacian: V Ат x (A, + Az + sr «7 4A.) (G q) 
7. 
Jacobian: JJ (5,8) as | (A As) (E, В.)-(Аг-А4) (8 - 8) (C-2) 
pe Aer B= YA: PB = 27, [(0-A9(8.-B)e (87 Ay) (Ba Ba) ез» 


The map factor "m" in the above equations is dependent on latitude 2 


and is expressed as: 


me SING, (с-4) 


|+ Sin 60 


The distance between gridpoints, "а" іп the grid used (Figs. 2, 


31), is 381 kilometers, true at 60 north. 
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Finite Differencing Grid 


Figure 31 
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APPENDIX D 


Computations of Static Stability Parameter 


The static stability parameter, Z , at level k may be expressed 
as follows: 


la as 


б = еті 


o 


which may be rewritten as: 


а= «ος μη) o 


by use of the Poisson relationship. Using the relation 


ΟΖ 


pr Ж. 30 (D-3) 
Ҡа Olu 
gives: 
ya 
= م‎ ‚_ Ra 90 (D-4) 
PS) Cw 22. 
Assuming constant within the layer pair (k-1, k), (k, k+1) 


the previous equation is set up as an ordinary differential equation in 
the form 
2. 
ο ὢ ἐν... а, =" 
T 6, Abe A ке 
d (In) e Abe 
The above is integrated with respect to Lap from Ber to Ze 


where these pressures are the log means for the two layers, that is: 


^". UE) m Au SR) 


3 L 
2 2 “ 


17 


After the above operation has been performed and solved for Oy, s it 


leads to the finite-difference form: 


4 m Фь Or Qu. + Dr. Фа , Eg (2,- 9. ) (D-6) 
24; 8,)) hk κ pum δι 
+! №. 


1 Ва (6 - 
2 2Cp СЖ $, 


The lower bound on the braced quantity above is to be 


thus (í. may be solved at all required levels. 
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APPENDIX E 


Procedure for 100-mb Height Extrapolation 


Navy Weather Research Facility (Lea, 1961) outlines a method of 
deriving a 100-mb extrapolated height field from the 200-mb heights 


and temperature fields in the form 


p: DX ANZ au ec (E-1) 


/οο 


The coefficients A and A. vary with both latitude and month of 


0” А» 2 


the year. 


0? А.» A, are defined for the centers 


of latitude belts of 10° latitude (or multiple thereof) in width. Since 


The regression coefficients A 


the A. for one belt generally is different from that of the adjacent 
belt, continuity required a smooth transition from one belt into the 
next, which was accomplished by linear interpolation between As and 
the linear distance between latitude band centers. 

Since temperature fields were not utilized directly in this study, 


the following approximation was made. 


2 $ 


| 
τ. yy (E-2) 
Τη τα ὃν 
Using the centered finite-difference scheme for o9 at 200 mb, 
Oln p 
equation (E-1) becomes 
8,22 
= — +273) 
= As EE Er Az Ra m3 (E-3) 


(| ع‎ 
Ra). 3 
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APPENDIX 


FORTRAN IV PROGRAM 
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